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' We provide a methodology for generating interatomic potentials for use in classical molecular 

, dynamics simulations of atomistic phenomena occurring at energy scales ranging from lattice vibra- 

■ tions to crystal defects to high energy collisions. A rigorous method to objectively determine the 

shape of an interatomic potential over all length scales is introduced by building upon a charged- 
ion generalization of the well-known Ziegler-Biersack-Littmark universal potential that provides the 
short- and fong-range hmiting behavior of the potential. At intermediate ranges the potential is 
smoothly adjusted by fitting to ab initio data. Our formalism provides a complete description of the 
I interatomic potentials that can be used at any energy scale, and thus, eliminates the inherent ambi- 

guity of splining different potentials generated to study different kinds of atomic materials behavior. 
We exemplify the method by developing rigid-ion potentials for uranium dioxide interactions under 
conditions ranging from thermodynamic equilibrium to very high atomic energy collisions relevant 
O . for fission events. 

Molecular Dynamics (MD) simulations provide a convenient tool for studying the dynamics of large atomic ensem- 
bles, provided that the dynamics of interest is not of a duration that makes the simulation time impractical. There 
are many examples of beautiful applications of how MD has been able to provide detailed insight into the dynamics 
and statistics of materials behavior. One contemporary interest is the field of high energy radiation damage of crys- 
talline materials, such as nuclear fuel. A core material of interest in this context is uranium dioxide, and one of the 
C important aspects of the interest in this material is to understand the evolution and statistics of atomic displacement 
, cascades due to high energy radiation^. Classical Molecular Dynamics is ideally suited for this kind of study since it 
' strikes a fine balance between being coarse enough to simulate the spatial scale necessary to represent the extent of 
O \ a damage cascade due to, e.g., a 100 MeV atomic collision with being detailed enough to retain the atomic structure 
of the material. The high energy range of the potential (short range in atomic separation) is consistent with the 
well-accepted Ziegler-Biersack-Littmark (ZBL) universal pair potentials^, which treat close range atomic interactions 
I as screened Coulomb forces between the nuclei. However, the complexity of the true interatomic interactions cannot 
^ ' be fully represented in an efficient manner by a simple classical functional form. Thus, one needs to develop a set 
OO , of essential interaction features that are necessary for a given application. This is a particularly challenging exercise 
for radiation damage simulations due to the disparate scales of energies involved. Therefore, interatomic potentials, 
suitable for this purpose, are typically constructed by smoothly joining different types of interactions. At medium to 
long range distances, a traditional potential (e.g., Buckingham, electrostatic etc.) fitted to a variety of thermodynamic 
and structural data is used. At short-ranges, accurate potentials are developed by fitting to ab initio data. The ZBL 
universal potential is one such very popular pair-potential developed by Ziegler et at in the 1980's, as a generic 
function of the atomic numbers of the species involvedi^ Although each type of interaction is directly determined 
through fitting, the determination of a suitable spline that smoothly joins these two pieces is a highly non-unique 
process. Splining leads to an inherent ambiguity in the behavior of the complete potential, since the exact cutoff 
distances and the spline's algebraic form will have consequences for how large the "cores" of the atomic interactions 
are and how the potential behaves in the region of transition. This, in turn, will have a significant impact on the 
ion trajectory and damage production one is ultimately interested in. Such ambiguity is illustrated in Figure [T] Of 
course, a complete description of the dynamics of the system will also require the inclusion of a suitable electronic 
stopping model^!^!^!^!^!^. 

In the present communication we introduce a rigorous method to unambiguously determine the form of the potential 
over all distances using ab initio data. Although our approach is applicable to any material system, we illustrate it 
using the important example of the nuclear fuel UO2, which has been the focus of several detailed computational 
investigations, both through ab initio^ and MD^ i^°'^^ methods, due to its critical importance in the nuclear industry. 

Our approach is to first generalize the universal ZBL potential to include charged ions that behaves correctly in 
both short range and long range limits. The advantage of building upon the ZBL formalism is that our potential 
automatically inherits the well-tested ability of the ZBL potential to describe high-energy scattering phenomena 
associated with the short-range behavior of the potential. This generalized potential smoothly interpolates between 
these two regimes over a physically-motivated length scale that is based on atomic orbital sizes instead of necessitating 
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a user-specified transition radius for the electrostatic interactions at short ranges to prevent double-counting. The 
only component that remains to be determined by fitting is then the medium-range energy contribution associated 
with chemical bonding. As this contribution is only significant over a relatively small range of distances it is possible 
to introduce lower and upper cutoffs, where this contribution must smoothly vanish. Importantly, these physically 
motivated cutoffs can be fitted to ab initio data and are therefore no longer arbitrary; unlike in the current practice. 

The ZBL potential^, which properly accounts for the screening of nuclear charge by the electronic clouds as a function 
of interatomic distances, is built on considering two interacting spherically symmetric rigid electron clouds. In this 
spirit, we also consider two interacting spherically symmetric rigid electron clouds with electron densities determined 
ab initio and fitted to a sum of Slater functions This approximation is valid beyond distances where electron clouds 
overlap and chemical bonds form. For short distances (i.e., distances less than chemical bond lengths) we obtain the 
energy as a function of distance between any two atoms through first order perturbation theory. It was demonstrated 
by Ziegler et al. that more sophisticated self-consistent field calculations incorporating the distortion of electronic 
clouds did not lead to any significant differences in the resulting interatomic potential at short distancesi^ Thus, since 
we employ the same electronic density and the same screening function (ratio of the actual atomic potential at some 
radius to the potential caused by an unscreened nucleus) as used by Ziegler et at, we recover the ZBL potential at 
these short distances, as we will show later. 

We consider two spherically symmetric charge densities per unit volume pi{r) and P2{i'), with central point charges 
of Zie and respectively, pi{r) and p2{f') being normalized to equal {Zi -I- qi)e and {Z2 -\- q2)& respectively. We 
further think of the point charge Zie as being made of two point charges, {Zi -\- qi)e and (— gi)e; is similarly 
decomposed into {Z2 + q2)e and {—q2)e. qie and q2e here denote the net ionic charges on atoms 1 and 2 respectively. 
We make this decomposition so that the Coulombic interaction term naturally arises in the expression for the net 
interaction potential V{r), 

2 

V{r) = ZBLz,+,,^z,+,, (r) + |^ +t,+t2, (1) 

where Z BLzi+qi.Z2+q2i''') denotes the ZBL form of interaction between two neutral atoms having atomic numbers 
Zi + qi and Z2 + q2, but using the screening length for Zi and Z2 ■ This is because the screening length is governed by 
the electronic charge distribution close to the nucleus and not far away from it. Since the extra charges qi and (72 have 
been added to the valence shells of neutral atoms with atomic numbers Zi and Z2, we do not change the screening 
length. This was found critical in recovering the standard short-range ZBL potential, ti denotes the interaction 
between the point charge —qie and the {Z2 + q2)e point charge plus electron cloud system of the atom 2, and is given 
by 

h (r) = - - /' 4^*'/'2 is)ds- r Ansp2 is)ds] , (2) 

while t2 is the converse remaining point-atom interaction, expressed similarly. Eq.(l) is correct for atomic separations 
smaller than the 'chemical bond length' (where it recovers the original neutral atom ZBL) and for very large atomic 
separations as well. In the latter case, only the 2nd term in Eq.(l) survives as we show below. 

The task now is the determination of poir) and pu{t)- It is here important to notice that the potential for very small 
interatomic separations is only as good as the ZBL form (see Eq. (1)). Thus, we use the charge densities employed for 
ZBL, which are primarily Hartree-Fock-Slater atomic distributions for most of the atomic pairs. We fit the numerical 
data for charge density used by Ziegler et al. to a sum of Slater functions. While the density p{r) is known to be a 
monotonic decreasing function of radial distance for all atoms, the graph of Attt^ p{r) exhibits a number of peaks (see 
Figured]) corresponding to atomic shells. To ensure the best possible accuracy, wc fit to Aur^ p{r) because this is the 
quantity entering Eq.(l). We find that 2 and 4 Slater functions are sufficient to capture the behavior of Anr"^ p{r) for 
neutral Oxygen and Uranium atoms (see Figure [2]), i.e., 

po(r) =aie-'=i'- + a2re-'=^'' (3) 



pu{r) = 6ie^'i'' + b2re-^^'' + bsr^e'^^'' + biV^e'^^'' . (4) 

Fitting the ZBL charge density with Slater functions does not exactly ensure that the areas under the 4Trr^p{r) curves 
for Oxygen and Uranium are 8 and 92 respectively. The Slater function fits therefore needed a slight modification in 
their pre-factors. In addition, the pre-factors in the above two equations as reported in Table U have been multiplied 
by 10/8 for Oxygen and 88/92 for Uranium because we are interested in the electronic cloud of the ionic species and 
not of the neutral atoms. We also experimented with more sophisticated corrections (e.g., using noble gas densities 
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instead) but this did not change the results by more than the intrinsic accuracy of the ZBL potential. By using Eqs. 
(3) and (4) in Eq.(l) and performing the integrations, we obtain the following pair potentials for Oxygen- Oxygen, 
Uranium-Uranium and Oxygen-Uranium respectively: 

Voo{r) = ZBLio^ioir) + foo{r)\ (5) 

47160?" 47reo r e 

Vuuir) ^ ZBL,,Mr) + Mll^! + _ ^lf,,^r)] (6) 

'Ineor 47reo r e 

,r , ^ ryr^T (4)(-2)e2 26^88 47r^ ... 46^ 10 4^^ ... 

Vouir) = ZBLssM^) + ^ 1 fuuir)] + -. [ foo{r)] (7) 



where we have 



= ^ - ^^[120 + m,r + Se^lr^ + 81^ + iy] + 
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2463 he-^^\ , , 



^ - + 4/.r + -f ^[2 - 2e-'- - /ire"'-] (9) 

We illustrate in Figure [3] how closely the potentials given in Eqs. (5) - (7) match, for small r, the neutral atom 
ZBL, and for large r, the relevant Coulombic interaction. 

As for any empirical potential, there is an intermediate distance range for which the interactions follow neither 
a ZBL nor a purely Coulombic form. A correction term is thus needed for this regime. We find this correction 
term by fitting to an extensive database of GGA+U ah initio calculations on UO2. GGA+[/ is known to provide 
electronic and magnetic behaviors of UO2 that are consistent with experiments':- ^ and a correct treatment of the 
localized and strongly correlated 5/ electrons of Uranium.^"* Our ah initio calculations also take into account the 
experimentally observed noncollinear antiferromagnetic magnetic moment ordering and the Oxygen cage distortion in 
U02'^ Therefore, in addition to capturing correct elastic and defect properties, our potential also covers a much more 
vast energy landscape in the material due to the richness of the ah initio data used. We fit to a database obtained 
by GGA-\-U calculations with the projector-augmented-wave method implemented in the VASP— package. In the 
GGA-I-J7 approximation, the spin-polarized GGA potential is supplemented by a Hubbard-like term to account for 
the strongly correlated 5/ orbitalsJ^ We use the rotationally invariant approach to GGA-I-J7 due to Dudarev et al^, 
wherein the parameter U-J is set to 3.99 eV. This is the generally accepted value for this parameter to reproduce 
the correct band structure for U02.''^ The magnetic moments were allowed to be fully non-coUinear. The ah initio 
database so obtained comprises: 

(i) isochoric relaxed runs on a 12 atom unit cell which was isometrically contracted and expanded by various 
amounts (i.e., equation of state calculations wherein each data point was calculated under constraint of constant cell 
volume), and for which an energy cutoff of 500 eV and a 8x8x8 k-point grid were taken; k-point convergence was 
ascertained before choosing this value for the k-point grid. The cell was allowed to relax in shape but not in size. 
Ionic relaxations were carried out until residual forces less than 0.01 eV/Awere achieved. 

(ii) static (i.e., no ionic relaxation) runs on 96 atom 2x2x2 supercell in which one atom at a time (i.e.. Oxygen 
or Uranium) was perturbed from its equilibrium position by varying distances in different directions. Energy cutoff 
was 500 eV. After performing convergence studies on the k-point grid, gamma point only version of VASP was found 
to be satisfactorily accurate for this. Note that any interactions between atoms and their periodic images do not 
systematically bias the fit of the potential because the same supercell geometry is used in both the ah initio and the 
empirical potential energy calculations. 
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(iii) 96 atom 2x2x2 supercell for the formation energies of three kinds of stoichiometric defects, namely Oxygen 
Frenkel pair, Uranium Frenkel pair and Schottky trio. The vacancies and the interstitials were taken as far from 
each other as the superceh would allow. The details of the calculations are the same as that for case (ii) above. 
Correct prediction of defect energies has been given great importance in generating interatomic potentials for cascade 
simulations. 

(iv) first-order transition states in a 2 x 2 x 2 supercell for the migration energy of Oxygen and Uranium vacancy and 
interstitial. Nudged Elastic Band(NEB) method^" in conjunction with the climbing image method^^. for determination 
of saddle point energy, as implemented in VASP, was used for this. 

With the ab initio database so generated, we now fit the final potential forms as follows. 

f4)f4)e^ 8e^ 88 in 
Vuuir) = ZBLssMr) + ^-^p^ + - -fuu{r)] (10) 



ZBL.oMr) ^[f - f /ooM] < r < 1.17A 

(-2)(-2)e^ J 5*^ order polynomial 1.17A < r < 2.28A 

~ ineor ] 3'''^ order polynomial 2.28A < r < 2.84A 

-603.268eVAVr^ r > 2.84A 



< r < 1.42A 
1.42A < r < 1.70A 
r > 1.70A 

The long-range Coulomb terms in the interaction here were calculated through the standard Ewald summation 
technique. The upper cut-offs for all terms except the Coulombic in Eq. (10) may be chosen as per the availability 
of computational resources. As seen from Eq. (10) we now have absolutely no splines for the U+'^-U^^ interaction, 
refiecting the fact that no chemical bonding takes place. There are splines in the other two interactions but these are 
now unambiguously determined since the respective cut-offs are not imposed but instead determined through fitting. 
The splines maintain continuity through the second derivatives of the potential and the specific form of the splines 
can be uniquely recovered from these conditions. Since they have been fitted to accurate ab initio data, these splines 
do not introduce any spurious wriggles in the potential. For the interaction between two Oxygen ions, the potential 
has one (and only one) minimum at rmm = 2.28A, as may be seen from Figure [H We have thus minimized the 
unphysical features of all of the interatomic potentials in UO2 , which were demonstrated for the particular case of 
Oxygen-Oxygen interaction in Figure [T] 

The downhill simplex method of Nelder-Mead was used to carry out the fittingi^ The fitting involved minimizing 
the sum of the squares of the differences between the ab initio energies and the energies predicted by the potential for 
all the classes of data points as detailed above. The package GULPS§ was used for energy calculations and for atomic 
positions optimization. The quality of fit for the equation of state data and the perturbed atom data can be seen in 
Figure [5l The ab mi<Jo/experimental and predicted defect formation/migration energies are compared in Table [Til 
while Table IIIII lists the predicted ground state lattice parameter and other elastic properties as compared with the 
corresponding experimental valuesi^i^^ (extrapolated accordingly) and with values obtained with the Morelon et al. 
potential^. The agreement is very satisfactory. 

As a final validation of the developed potential we considered various dynamic properties by performing MD 
simulations in a constant number, pressure and temperature (NPT) ensemble comprising 6x6x6 unit cells. The 
system was equilibrated for 10.0 ps, while production runs were carried out for 100.0 ps, with a time step of 0.001 
ps and sampling every 0.05 ps. The fiuorite structure remained stable during all the runs we performed, up to 
temperatures of 2500 K. The properties we considered are the variation in the lattice parameter and the enthalpy 
as functions of the temperature. These are compared in Figure [HI with the corresponding experimental datai^ The 
quality is similar to what is given by the previous potentials22ii^'22i, as tabulated in the work by Morelon et al^ 

To summarize, we have shown a methodology for developing an interatomic pair potential such that it is appropriate 
for all relevant interatomic separations, without the need for any ambiguous splines. Splining between regions of 
different characteristics is not just an inconvenience in the implementation of a potential in MD simulations, but also 
introduces an uncertainty regarding which distances, and by which functions, one realizes the splines. The potential 
we have focused on in this presentation, the nuclear fuel UO2, is just one example of a model material in which 
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very relevant materials physics depends on accurate and reliable interactions over many orders of magnitude, and 
we have obtained the first complete description that allows for direct simulations of damage cascades due to high 
energy radiation effects. The potential has been generated based on a slight revision of the ZBL universal potential 
to account for ionic materials with the intermediate interatomic distances fitted to a broad database of ab initio 
structural energies. In view of these qualities, we expect it to be a very reliable potential for studying displacement 
cascades in UO2. 
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TABLE I: Values of coefRcients in Slater functions in Eqs. (3) and (4) 

ai = 2799.625 eA"^ &2 = 211.038 eA"* ki = 30.76 A~^ 

k2 = 6.77 A"^ bi = 3092188.94 eA'^' b2 = 13255095.09 eA"'' 

ba = 4982192.00 eA''^ hi = 135624.70 eA"'' li = 309.92 A'^ 

I2 = 87.23 A"^ I3 = 32.98 A'^ U = 13.80 A"^ 



TABLE II: Comparison of defect formation and migration energies (all in eV), between our values and best values as per 
previous potentials;^'' compared with ah initio values from this work and with experimental values.— 





Exptl.(E)/ ah miUo(Al) 


This work 


Previous Potentials 


Oxygen Frenkel Pair Formation Energy 


3.5 +/- 0.5(E), 3.9 (AI) 


3.3 


3.17 


Uranium Frenkel Pair Formation Energy 


9.5 -12(E), 10.1 (AI) 


15.5 


12.6 


Schottky Trio Formation Energy 


6.5+/- 0.5(E),7.4 (AI) 


7.1 


6.68 


Oxygen Interstitial Migration Energy 


0.9-1.3(E),1.5(AI) 


1.4 


0.65 


Oxygen Vacancy Migration Energy 


0.5(E),0.9(AI) 


0.5 


0.33 


Uranium Interstitial Migration Energy 


2.0(E), 1.3(AI) 


2.3 


5.0 


Uranium Vacancy Migration Energy 


2.5(E),2.8(AI) 


2.8 


4.5 



TABLE III: Comparison of various ground state elastic properties, between our values and best values as per previous 
potentials,— compared with (extrapolated) experimental values.^ N.B. These are predicted and not fitted values. 





Exptl. 


This work 


Previous potentials 


Lattice Parameter (A) 


5.46 


5.46 


5.46 


Bulk Modulus (GPa) 


207 


210 


125 


Elastic Constant Cn (GPa) 


389.3 


401.8 


216.9 


Elastic Constant C12 (GPa) 


118.7 


114.1 


79.1 


Elastic Constant C44 (GPa) 


59.7 


107.8 


78.5 
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1.2 1.4 1.6 1.8 2 . 2.2 2.4 2.6 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 

distance between atoms(A) distance between atoms (A) 




distance between atoms (A) 



FIG. 1: (a) ZBL screened potential (solid line) and Morelon et al. potential (dashes), joined together by a 5th order polynomial 
(plus signs), for the case of two Oxygen atoms, (b) First derivative of the net potential resulting from the same, (c) Second 
derivative of the potential. Since the spline was not fit to any data, one can not decide whether the resulting behavior in (b) 
and (c) is correct or spurious. 




FIG. 2: Charge density times 47rr^ fitted to sum of Slater functions for (a) Oxygen and (b) Uranium. Open circles denote 
values used by Ziegler et al.,^ while the solid lines indicate our fit using sum of Slater functions. The resultant error in the 
short range interatomic potential as compared to ZBL's original potential was well within the latter's standard deviation for 
both (a) and (b). Thus trying to capture more peaks for Uranium, by introducing more Slater functions, was not necessary. 
The coefficients used here are provided in Table HI 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.5 2 2.5 3 3.5 4 4.5 5 
distance between atoms (A) distance between atoms (A) 



FIG. 3: Comparison of our analytical potential form (dashed line) with (a) the currently used neutral atom ZBL interaction 
(solid line) for small distances, and (b) the ionic coulombic interaction (solid line) between two (-2e)point charges for large 
distances for the case of two Oxygen atoms. 




1.2 1.4 1.6 1.8 ^ 2 2.2 2.4 2.6 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 

distance between atoms (A) distance between atoms (A) 




1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 
distance between atoms (A) 



FIG. 4: (a) The interatomic potential for Oxygen-Oxygen as per current work (Eq. 10). (b) First derivative of the net potential 
resulting from the same, (c) Second derivative of the potential. These are to be compared with Figure 1. 
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FIG. 5: Quality of fit from our fitted potential for various ab initio energies: (a) expansion/contraction (open circles), (b) 
Oxygen atom perturbation (plus signs), (c) Uranium atom perturbation (cross signs). Asterisks denote experimental data. For 
each of Oxygen and Uranium, the first four perturbations are along < 100 > direction while the second four are along < 110 > 
direction. 




Temperature (K) Temperature(K) 



FIG. 6: (a) Relative lattice parameter variation using potential from current work (open circles) compared with corresponding 
experimental values (dashed line).^^ The scatter in the experimental values is also shown (solid lines), (b) Enthalpy variation 
using potential from current work (open circles) compared with corresponding experimental values (solid line).— The scatter 
in the experimental values was less than 1 percent and is thus not shown. 



